Force correlations in the q— model for general q— distributions 
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We study force correlations in the g-model for granular media at infinite depth, for general q- 
distributions. We show that there are no 2-point force correlations as long as q-values at different 
sites are uncorrelated. However, higher order correlations can persist, and if they do, they only 
decay with a power of the distance. Furthermore, we find the entire set of g-distributions for which 
the force distribution factorizes. It includes distributions ranging from infinitely sharp to almost 
critical. Finally, we show that 2-point force correlations do appear whenever there are correlations 
between g-values at different sites in a layer; various cases are evaluated explicitly. 

PACS numbers: 02.50.Ey, 45.70.Cc, 81.05.Rm 
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I. INTRODUCTION 

One of the main challenges of granular media is to char- 
acterize the network of microscopic forces in a static bead 
pack. In order to describe the corresponding force fluc- 
tuations, Liu et al. Q introduced the g-model. In this 
model, the beads are placed on a regular lattice and the 
(scalar) forces are stochastically transmitted, by random 
fractions denoted by the symbol q. Even in its simplest 
version, where one assumes a uniform q-distribution, it 
already reproduces the main feature of the experimen- 
tal observations: the probability for large forces decays 
exponentially jj], ||, |!|. Although for this uniform q- 
distribution the forces become totally uncorrelated, in 
general, correlations do persist j|. In the present study, 
we investigate for which ^-distributions this is the case 
and we reveal the surprising nature of these correlations. 
In order to perform an analytical study, we restrict our- 
selves to the scalar g-model and allow only correlations 
between g-values in a layer. More sophisticated lattice 
models, which include the vector nature of the force and 
allow correlations between layers are not considered here 

i- 

Although the q model is particularly simple, its be- 
havior turns out to be very rich. First of all, there is a 
so-called critical q-distribution, which produces a force 
distribution that decays algebraically instead of exponen- 
tially j|, ||. It therefore forms a critical point in the 
space of g-distributions, and its properties were recently 
investigated in great detail A second intriguing 

issue concerns the top-down dynamics of force correla- 
tions (the downward direction can be interpreted as time) 
J?], ||, H . Even if both in the initial state (top layer) and 
in the asymptotic state (infinite depth) all forces are un- 
correlated, there will be correlations at all intermediate 
levels. Correlations become longer in range while their 
amplitudes diminish in a diffusion process, and as a re- 
sult, the asymptotic force distribution is only approached 
algebraically ||. This process is closely related to the 
subject of this study, namely the presence of force corre- 
lations at infinite depth. 

Let us recapitulate the definition of the g-model. The 
beads are assumed to be positioned on a regular lattice. 




FIG. 1: The displacement vectors a in the q-model for (a) 
the triangular packing (side view) and (b) the fee packing (top 



Let fi be the force in the downward direction on the ith 
bead in a layer. This bead makes contact with a num- 
ber of z beads in the layer below, which we indicate by 
the indices i + a. The cv's are displacement vectors in 
the lower layer as shown in Fig. [l]. Bead i transmits 
a fraction qi^ a of the force fi to the bead i + a under- 
neath it. These fractions are taken stochastically from a 
distribution satisfying the constraint 



(1) 



which assures mechanical equilibrium in the vertical di- 
rection. So, we can write the force /' on the jth bead in 
a layer as 



fj — Qj-a,a fj 



(2) 



As the weights of the particles are unimportant at infinite 
depth, we have left out the so-called injection term. The 
distribution of forces at infinite depth depends on the 
q-distribution H(q), where the symbol q is a shorthand 
for all the <7i )Q at a given layer. This H{q) can be any 
function that is constrained by Eq. (|l|) . If we now assume 
that there are no correlations between the g-values at 
different sites, the q-distribution is of the form 



i \ a / 



Qi = {Qz.a}, (3) 



where T](qi) is symmetric in its arguments qi, a - Although 
we will refer to these (/-distributions as "uncorrelated" , 
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note that there are always correlations between the qi^ a 
of the same site due to the 5-constraint. 

In the first part of this study, we show that there is 
only a limited set of T)(qi) for which the stationary force 
distribution can be written as a product of single-site dis- 
tributions, and therefore is totally uncorrelated. This set 
is an extension of the set that was already identified by 
Coppersmith et al. In their extensive study, they also 
provided numerical evidence that, in general, correlations 
can persist. We will show that correlations are still ab- 
sent in the second order moments. However, higher order 
correlations do exist and surprisingly enough, these turn 
out to decay algebraically. The results for the triangular 
packing and the fee packing are summarized in Table |, 
section VII . In the last part of this work, we show that 
one induces 2-point force correlations by allowing cor- 
relations between (/-values on different sites in a layer. 
These correlations will generically vanish with a power 
law, except for the triangular packing, where the decay of 
force correlations follows the decay of the g-correlations. 

The paper is organized as follows. In section |l| we 
derive a criterion that a distribution rj{qi) has to obey 
in order to produce an uncorrelated stationary state. We 



then show in section III , that this criterion is only obeyed 
for a limited set of rj(qi). After that, we study the na- 
ture of the correlations, by writing the evolution of the 
force moments as master equations in section [V , and by 
analysing the stationary solutions of these equations in 
section M. Section VI deals with the effects of allowing 
correlations between the ql of different sites in a layer, 
and the paper closes with a discussion. 



II. CRITERION FOR FACTORIZATION 

Using the recursive nature of the force transmission, 
Eq.(^|), one can write down the following recursive rela- 
tion for the force distribution 0, 01 : 



P'(f') = j H(q)dqJ P(f)df 

x n - n fj-°)> w 

j a 

where we have introduced a vector notation for the forces 
in one layer / = (/i, • • • , /jv)> and for the integrations we 
use the abbreviations 



yv=nf * 



(5) 



/ d< i= n / d * = n n / d ^°- ^ 

J i J i a J o 

It is often convenient to work with the Laplace transform 
of Eq.(Q). Defining the Laplace transform as 



the recursion simplifies to ^, |TJ] 

P'(s) = f H{q)dqP{s{q)), 



with 



Si{q) = y^gi.q Si + c 



(8) 



(9) 



The two representations Eq.([|) and (ra) are equivalent, 
and they will both be used, depending on the nature of 
the problem. 

The force distribution at infinite depth P*(f) or P*(s) 
can be obtained by finding the fixed point of the recur- 
sive relation. The main question of this section is to 
determine whether a given H(q) leads to a P*(f) that is 
simply a product of single-site force distributions p*(fi). 
In section VI we will show that this can only be the case 
for q-distributions of the type Eq.(||). So for this sec- 
tion, the question is: which rj(ql) lead to uncorrelated 
asymptotic states? 

To answer this question, let us assume that such a fixed 
point exists, i.e. 

p*(j) = n p*(/i)> ° r p *w = n ( io ) 



Inserting this Ansatz into the Laplace representation of 
the recursion relation, Eq.(B), yields 



P*(s) 



II / rM)6 \1 -5>,„j dq t p* 



(11) 



where the function ip(si+ ai , ■ ■ ■ , Si+ az ) is the outcome 
of the integral over the ql. The arguments represent 
the z sites that are connected to site i in the previous 
layer. Integrating out all forces except those at the z 
sites connected to i means putting all Sj = except the 
set {s i+a }: 

P (^i+ai i ' ' ' ? Si-\-a. z ) — ^(^i+CKi •> ' ' ' j <~>i-\-a z ) 

x Y[i>( Si+a ,o,---y- 1 .(i2) 

a 

This projection of the total force distribution can only 
factorize if ip(si+ ai , • • • , Sj+a z ) is a product function as 
well, i.e. 



(13) 



P(s)= / dfexp(-s-f)P(f), 



(7) 



This leads to the following criterion for asymptotic fac- 
torization: 
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• Given a (/-distribution 77(g), one can construct a 
factorized fixed point if, and only if, there is a func- 
tion ip(s) that satisfies the following condition: 



dq 



a 

(14) 



• This function ip(s) is related to the single-site dis- 
tribution as 



p*(s) = $( s ) 



(15) 



Here, we omitted the site index i, and furthermore, our 
formulation depends only on z (the number of g-values 
per site) and not on the details of the lattice. 



III. SPECIAL CLASS OF ^-DISTRIBUTIONS 
LEADING TO FACTORIZATION 

It is a well-known fact that the so-called uniform dis- 
tribution, in which rj(qi) is a constant, produces an un- 
corrected asymptotic force distribution. In fact, Copper- 
smith et al. identified a countable set of q— distributions, 
of which the uniform distribution is a member, that have 
this property B. Although it might seem obvious that 
a uniform distribution leads to an uncorrelated asymp- 
totic state, it is really not trivial. Due to the constraint 
of Eq.(Q), there are correlations between the qi <a on each 
site i, which induce force correlations that only disappear 
under the special conditions discussed in the previous sec- 
tion, Eq.(|14|). In this section, we will show when these 
special conditions are obeyed. 

There is a mathematical relation that is extremely im- 
portant for the (/-model [fol: 



(16) 



It holds for any real r > 0. From this relation, it is 
immediately clear that for all (/-distributions of the type 



T(zr) 



there is a ip{ s ) that obeys Eq.(|lJ), namely 

1 



n(?-r _i . r>° (17) 



(18) 



The corresponding single-site force distributions are 



p*(s) 



(1 + sj zr) z1 1 (zr) 

(19) 



We rescaled the Laplace variable s, in order to put 
(/) = 1. Coppersmith et al. already found these q- 
distributions for integer values of r, also based on Eq.(|l6|) 
[Q. However, it holds for any real r > 0. This means that 
the set for which the stationary force distribution factor- 
izes is substantially larger; it ranges from the infinitely 
sharp distribution (r — > 00) to the critical distribution 
(r — > 0) pi] . Note that one recovers the results for the 
uniform distribution by putting r = 1. 

Although there is a huge variety of (/-distributions that 
lead to uncorrelated force distributions, in general one 
cannot find a "4>(s) that obeys Eq.( |l4[) . We will prove 
this by making a Taylor expansion of ip( s ) 



(20) 



n=0 



and then try to solve for the coefficients ip n by impos- 
ing Eq.([l4j). It turns out that the equations can only 
be solved under special conditions, which are precisely 
obeyed by the class of (/-distributions given by Eq.(^). 

Let us first focus on the left hand side (LHS) of Eq.([i"4"|). 
The Taylor expansion will give rise to terms of the type 
(qisi) ni ((/2S2)™ 2 • • • (qzS z ) n *, which have to be integrated 
over all q a . This leads to terms s 7 ^- sVj 2 ...s™ z with prefac- 
tors given by the moments of r](q) 

q^qfZqF = J ^6 (l " I>1 dq <HHT ■■■<&* ■ 

These moments are not independent, due to the con- 
straint Eq.(|l|). In appendix |A|, we show that the mo- 
ments 



(22) 



are in fact sufficient to characterize all relevant moments 
of Eq.(pl|). Besides the moments, there are of course 
additional prefactors consisting of combinations of the 
ipn] these are the quantities we try to find, for a given 
(/-distribution rj(q). 

The right hand side (RHS) of Eq.Q) also produces 
terms s^s? 2 ■ ■ ■ s" z , with prefactors tp nx tp n2 • ■ -ipn^- The 
remaining task is to equate the prefactors of the terms 
s i ls 2 2 ' ' ' s " z on both sides of the equation. This gives a 
set of equations, from which one can try to solve for the 

Ipn- 

The zeroth order equation is trivially obeyed for any 
■00 j as can be seen by putting all s a = 0. For convenience 
we fix ipo = 1. The same happens at first order, since for 
each a, the LHS contains z terms ipiq^s a — l/zipxs a , 
and the RHS is simply ipi s a- The first non-trivial equa- 
tion appears at second order. There are two equations, 
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for si, and for s a s c 



Z1p2 



where a ^= a': 

Zi>2+ 2 - M)V2 

z(z-l) (2 \ 2{l-z m ) 



z(z-l) 



1p2 
V>1 



(23) 



Due to the constraint ^2 a qi, a = 1, one can obtain 
an identity by multiplying the first equation by z, and 
adding it to the second equation multiplied by z{z— 1)/2. 
Hence, the two equations are not independent and i/>2 can 
be solved. The value of ip2 depends only on 772, the second 
moment of the g-distribution [jjjl . 

Working out the combinatorics of the higher orders, 
one finds the following general mathematical structure: 

• At the nth order, there are as many equations as 
there are different partitions {ni, 77,2, • • • , n z } that 
make n a — n. Permutations should not be 
considered as different because 77(g) is symmetric 
in its arguments. 

• One of these equations is dependent, as one can 
obtain an identity by adding the equations, after 
multiplication by appropriate factors. 

For z = 2, there are two third order equations, corre- 
sponding to the partitions {3, 0} and {2, 1}, of which only 
one is independent. This means that ips can be solved as 
a function of 772 (in appendix |A| we show that 773 depends 
on r/2, for z — 2). We run into problems at fourth order, 
where we have {4,0}, {3,1} and {2,2}, and hence two 
a priori independent equations for one coefficient 1^4. It 
turns out that the remaining equations are only identical 
if there is a relation between 774 and 772, namely 



3077I - 11772 + 1 
16772 - 2 



(24) 



In appendix |A|, it is shown that this relation is precisely 
obeyed by the class of q- distributions Eq. ( |l7| ) for which 
if)(s) was already solved. 

The fact that the expansion of tp(s) = [p*(s)] 1 ' z only 
fails at fourth order implies that a mean field approxi- 
mation, in which one explicitly assumes a product state, 
does give the exact results up to the third moment of 
p*(f). This is precisely the reason why the mean field 
solution p ' (f) differs only marginally from the real so- 
lution. To be more precise, the deviation p ml (f) — P*(f) 
should change sign 4 times, since it does not affect all 
moments lower than (/ 4 ). A careful inspection of the nu- 
merical results in Q for a g-distribution in which q = 0.1 
or q = 0.9 shows that these small "wiggles" are indeed 
present. To magnify this effect, we show our simulation 
data in Fig. ||. 

For z — 3, the problems already appear at third or- 
der. Since we have {3,0,0}, {2,1,0} and {1,1,1}, we 
encounter two independent equations for ^3. Again, it 




/ 

FIG. 2: Numerical simulation of a g-distribution with q = 0.1 
or q = 0.9. The small deviation p m (/) — p* '(f) changes sign 
4 times. 



turns out that the equations can be solved if there is an 
additional relation between the g-moments: 



"3 



15771 



9r/ 2 



(25) 



For z > 3, there are two independent third or- 
der equations as well, originating from {3,0,0,0,---}, 
{2, 1, 0, 0, ■ • ■} and {1, 1, 1, 0, ■ ■ •}. This problem can al- 
ways be overcome by assuming a particular relation be- 
tween the moments 773 and 772 , corresponding to the spe- 
cial g-distributions of Eq.(|lj). Since at higher orders the 
number of equations per coefficient ip n becomes increas- 
ingly high, there will be no other g-distributions than 
those of Eq. (|l7|) that obey Eq.(|l4|), and thus have an 
uncorrelated force distribution. 



IV. EVOLUTION OF MOMENTS 

Now that we know that, in general, correlations do ex- 
ist in the stationary force distributions, it is interesting 
to study the nature of these correlations. In this section, 
we write the evolution of the moments as master equa- 
tions, along the lines of Ref. 0. With this formalism, 
we will, in the next section, analyze the correlations by 
finding the stationary states of these master equations. 

First, let us define the second moments of a distribu- 
tion as 

M 2 (k) = ( fif i+k ) - j dfftfi+k P(f). (26) 

We have reintroduced the site-index i, and k is a 
displacement-vector in a layer. As the system is trans- 
lationally invariant, these second moments depend only 
on the displacement k. The recursion for these moments 
is obtained by combining Eq.(|^) and Eq.(||) as 



M ^) = E(y h ^<% 

Qj,a Qj-\-k-\-a — a' .a' 

a,a f 

x M 2 (k + a- a'). 



(27) 
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Using the overline notation for the g-averages again, 
Eq.(p7|) becomes 

M 2( k ) = ^ gj,a qj+k+a-a'.a' M 2 (k + a - a'). (28) 

a, a' 

This relation is reveals from which points (in correlation- 
space) the moment M 2 (k) receives a contribution during 
a recursion step. However, it is in fact easier to consider 
the opposite relation, revealing how much a moment con- 
tributes to correlation space points during recursion. The 
"inverse" of Eq. ( p8| ) becomes 

M 2 (k) — > qi. a q l+ k, a ' M' 2 {k + a' - a) for all a, a'. 

This latter relation allows for a master equation type 
formulation, as we may write it in the form 



E 

7 



M' 2 (k)-M 2 {k) = 
W 7 {k - 7) M 2 (k - 7) - W-y{k) M 2 (k). 



The transition rates are defined as 

W 7 (k) = qi, a q i+ k,a', 
with 7 determined by the set a, a' as 
7 = a' — a. 



(30) 



(31) 



(32) 



In the current problem, where we consider second order 
moments, the transition rates are particularly simple. If 
k 0, the g-averages are independent, and will always 
give the value 1/z 2 (this only holds for g-distributions of 
the type Eq.©). If k = 0, one encounters second mo- 
ments of T)(q), as in Eq.(|T]). This leads to the following 
transition rates: 



k = W (0) = r l2 W 7 #o(0) = 
fc^O => W 1 {k) = \. 



1 ~ ZT) 2 

z(z-iy 



(33) 



So, the moments evolve in an anomalous diffusion pro- 
cess, with differing transition-rates at the origin. For a 
detailed discussion of the corresponding dynamics, see 
Ref. U . Note that this diffusion takes place in a d — 1 
dimensional space, as a, and therefore also 7, is a dis- 
placement in a layer. In the remainder of this paper we 
use the bold notation 7 whenever the displacement is 
really a vector. 

The advantage of this somewhat formal representation 
is that we can take it over to higher order moments with- 
out further ado. The generalization of the master equa- 
tion for the nth order moments M n (r) becomes: 

M' n {v) - M n {v) = 
W^r-i) M n (r - 7) - W_ 7 (r) M„(r), (34) 



with the position indices r = (fci, k 2 , • • • , k n —i), and the 
displacements 7 defined as 

7 = (ai — a, a 2 — a, ■ ■ ■ , a«— 1 — a). (35) 

The dimensionality of the diffusion process has now be- 
come (n— l)(d— 1). The transition rates can be calculated 

as 



W 7 (r) = qi, a qi+k! 



(36) 



Analogous to the second moments, these transition rates 
are all l/z n , as long as the indices of the position vector r 
are not equal to zero nor coincide. However, the differing 
rates make the problem complicated, because one has to 
deal with different transition rates at special points, lines, 
planes etc. in the space of diffusion. 

One can now study the correlations at infinite depth by 
finding stationary states of the master equation for the 
moments. As a first attempt to construct a stationary 
solution, i.e. M' n {v) — M„(r) = 0, one can try a detailed 
balance solution. Detailed balance means that there is 
no flow of "probability" from one point to another. In 
that case, all terms of the sum on the right hand side of 
Eq.(p4T) vanish individually, i.e. 



W^(r)M n (r) 



W 7 (r • 



7 )M„(r- 7 ) all r, 7 . (37) 



This condition can also be formulated in terms of ele- 
mentary loops, which are the smallest possible pathways 
from a point to itself. For all lattices in this study, these 
elementary loops are triangles, and we denote the three 
jump rates as (a, 6, c) or (a', 6', c') depending on the di- 
rection in which the loop is traversed. It is easily verified 
that the property 



(38) 



must be obeyed in all elementary loops in order to have 
a detailed balance solution. In the next section we show 
that correlations appear whenever the detailed balance 
conditions are not obeyed. 



V. HIGHER ORDER CORRELATIONS 

In this section, we study the nature of the correlations 
for ^-distributions of the type Eq.© that do not fall into 
the special class of Eq. (p"7|) . We first solve the stationary 
master equation for the second order moments, for which 
we already know that there are no correlations (section 
III). For the triangular packing [z = 2), correlations only 
show up at fourth order, and these fall off as l/r 5 . For 
z > 3, there are third order correlations that also decay 
with a power law; for the fee packing (z = 3) the decay 
is l/r 4 . Finally, we provide a simple relation to calculate 
the various exponents. 
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A. Second order moments: no correlations 

In order to get familiar with the structure of the mas- 
ter equations, we first consider the second order moments 
desribed by Eq.(|3(]). Away from the origin k = 0, all 
transition rates of Eq.(|33"|) are identical. Therefore, the 
detailed balance condition Eq.(|37|) requires all M2(k =/= 0) 
to be identical. The value at the origin M2(0) has to 
obey a detailed balance condition for each 7^0, but 
these equations are identical for all 7 because the corre- 
sponding rates are the same. Putting M%{k 7^ 0) = 1, 
one obtains the following stationary solution: 



(/</<+*> = < *(!-*»») ' ■ (39) 

[1 , fc^o 

This solution precisely describes an asymptotic state 
without any 2-point correlations, as the average of the 
product (fifj) equals the product of the averages for all 
i 7^ j. Of course, any multiple of Eq.(^S|), also forms 
a stationary solution of Eq.(p0[). However, these so- 
lutions are physically irrelevant in the thermodynamic 
limit, where the lattice size — > 00 ||. Moreover, we find 
that the asymptotic second force moment is solely deter- 
mined by z and 772- For critical ^-distributions one has 
r]2 = 1/z, leading to a diverging second moment. 

B. Third order moments 

The diffusion of third order moments (fifi+kfi+i) takes 
place on a 2{d — l)-dimensional lattice, since there are 
two free parameters k and I of dimension d — 1. On this 
lattice, there are three special subspaces, namely k = 0, 
I = 0, and k — I, for which the transition rates of Eq.(|3^) 
differ from the bulk-value l/z 3 . Moreover, the rates at 
the origin k = I = differ from both the bulk-rates and 
the rates on the special subspaces. 

Let us first consider the triangular packing (z = 2), 
for which the third order moments diffuse on a 2- 
dimensional lattice, with differing rates on three special 
lines. As these lines are all equivalent, it is natural to 
draw them at an angle of 120°, see Fig. ||. We then 
obtain a triangular lattice, with transitions to six near- 
est neighbors and two self jumps, which are "transitions" 
to the same lattice site (7 = 0). The detailed balance 
condition between a special line and the bulk is natu- 
rally identical to the second order condition, implying 
the same ratio as in Eq.(p9|). As the transition rates at 
the origin are again identical for each 7 7^ (because 
of symmetry), one can construct the following detailed 
balance solution: 

{m . . 
TTW ' ° ngm 
I lines • ( 4 °) 
2(1-2^) ' 
1 , bulk 



This means that there are also no 3-point correlations for 
z = 2: at the origin we encounter (/ 3 ), on the lines we 
have {fffi+k) = (f 2 )(f), and in the bulk (fdj +k fi+i) = 
(J) 3 . It is easily checked that condition Eq. (pq) is indeed 
satisfied in every elementary loop. 




1=0 



FIG. 3: Triangular packing: third order moments diffuse on 
a triangular lattice. 

For the fee packing (z — 3), the third order moments 
diffuse on a 4-dimensional lattice. Unlike the z = 2 pack- 
ing, it is not possible to construct a detailed balance so- 
lution in this case. First, we write the displacement vec- 
tors as 7 = (a' — a, a" — a) = (71, 72), where the a's and 
7's are 2-dimensional vectors (Fig. [l]). One can jump 
away from the origin with two different rates, namely 
qfq2 and gTgiW- These rates correspond to 71 = 72 
(towards a special plane) and 71 7^ 72 (into the bulk) 
respectively. Checking the detailed balance condition in 
the elementarytriangle origin-plane-bulk-origin, it turns 
out that Eq. ( p8|) is only obeyed if 773 and 772 are related 
as in Eq. (p5|) . Of course, this is precisely the case for 
the class of Eq. ([l7]) for which we know that asymptotic 
factorization occurs. In general, however, it is not possi- 
ble to construct a detailed balance solution for the third 
order moments. In the next paragraph, we show that 
the absence of detailed balance indicates that there are 
force-correlations, that decay with a power law; in this 
case the decay is 1/r 4 . 

C. Fourth order moments 

The fourth order moments (fcfi+kfi+ifi+ m ) of the tri- 
angular packing diffuse on the bee lattice depicted in 
Fig. Q The three directions k,l,m precisely define a 
bec primitive cell [^3|. There are now differing rates at 
the origin as well as on lines and planes for which one or 
more indices coincide or are equal to zero. The precise 
geometrical structure is explained in appendix |^. There 
are now two a priori different directions away from the 
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origin, that is to corners (fffi+i) and to b ody centers 
(fffi+i)- Checking the loop condition Eq.([38|) for the 
loop origin-corner-body center-origin, one finds that it 
is only satisfied when 774 and 772 are related as in Eq.(24). 




corner 
body-center 



FIG. 4: Triangular packing: fourth order moments diffuse on 
a bcc lattice. 

The question that emerges is: What are the stationary 
solutions of the master equation, when the detailed bal- 
ance condition is frustrated at the origin? To answer this 
question we first consider a simplified version of the bcc- 
problem, as a first order approximation. In this simple 
version, we assume that all jump rates are 1/z 4 = 1/16, 
except at the origin where we distinguish between the two 
different directions. Although we neglect the differing 
rates on the special lines and planes, the loop condition 
is still frustrated in the elementary loop origin-corner- 
body center-origin. Using X^ 7 W 7 (r) = 1, we write the 
stationary master equation as 



M(r) = ^W r (r- 7 )M(r- 7 ) ) 



(41) 



[1 - 2Wo(r)] M(r) = £ W y (r - 7 )M(r - 7). (42) 

This allows us to eliminate the two self rates Wo by means 
of a simple transformation: 



M(r) = [l-2Wo(r)]M(r), 
W 7 (r) = W 7 (r)/[l-2W (r)] 



(43) 



The sum over the new rates again adds up to unity and 
Eq.([42"|) becomes 



M(r) = ^W 7 (r-7)M(r- 7 ). 

7#0 



(44) 



Hence we can omit the self jumps by first solving the 
equation for the "hatted" variables, and then transform- 
ing back to M(r). As Af(r) — > 1 for large r, it is conve- 
nient to write 



M(r) = - [l 



- SM(r) 



(45) 



The quantity 8M (r) is in fact the appropriate measure 
for correlations |1J] . After eliminating the two self rates, 
all jump rates have become 1/14, except at the origin 
where the rates to the 8 corners (c) can differ from the 
rates to the 6 body centers (b). We therefore have 



W-y{r) = 1/14 + <5(r)e. 



7- 



(46) 



The rates to the corners are denoted by e c and those 
to the body centers by £&. They fulfill the condition 
8e c + 6eb = 0. This results in the following equation: 



5M(r) 



1 

14 



^<5M(r 

7#0 



7) = -M(0)5> 7 <Kr- 7 ) 



7#0 



(47) 

Note that this is a discrete version of Poisson's equation: 
the LHS is a discrete Laplacian and the RHS, originating 
from deviating rates, acts as a multipole around the ori- 
gin. This equation is solved in appendix [b] by a Fourier 
transformation, leading to 



7 



M(r) = 1 + M(0) ^ i E{ ^ {k) expHk ■ r). (48) 

The functions D(k) and E(k) are defined in appendix 
[b|; 1 — D(k) comes from the discrete Laplacian (in the 
continuum equation it would simply be k 2 ), E(k) is the 
Fourier transform of the source, and the sum over k is 
the inverse Fourier transformation running over the Bril- 
louin Zone. The amplitude of the source M(0) can be ob- 
tained self-consistently, by setting r = 0. This involves a 
complicated integral over the Brillouin Zone (BZ) of the 
bcc-lattice; the outcome, however, will be of the order 
unity. The large r behavior of the correlations is deter- 
mined by the small k behavior, so E(k)/(1 — D(k)) has 
to be expanded around k = 0. The first term that gives 
a contribution is 



49e c f dk (k%k. 



24 



2J.2 



kjkx) exp(-ik • r) 



Vbz 
343e r 



32tt 



x 2 y 2 



- y 2 z 2 + z 2 x 2 



(49) 



The solution of Eq.(|47|) decays as 1/r 5 ; the terms x 2 y 2 
etc. give the proper angular dependence. This result 
can be directly understood from the analogy with elec- 
trostatics. The solution of Poisson's equation Eq.(|47|) 
can be expanded in asymptotically vanishing spherical 



harmonics: Y/ m / 



J+i 



The symmetry of the bcc lattice 



allows only harmonics with 1 
1/r 5 decay. 



> 4, leading to the observed 
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So we find that the stationary master equation for the 
moments becomes a discrete Poisson's equation, and the 
presence of differing transition rates leads to a multipole 
source around the location of these rates, see Eq.(47). 
However, this source is only "active" if there is no de- 
tailed balance, since detailed balance leads to trivial so- 
lutions like Eq. ([h]) . Keeping this in mind, let us now 
investigate the real fourth order problem, including the 
differing rates at the special lines and planes. We argue 
that the asymptotic value is still approached as 1/r 5 , but 
the amplitude of this field will be modified. Since there 
is no detailed balance, the differing rates at the lines and 
planes will act as sources as well. Their amplitudes, how- 
ever, will decay with increasing distance, since the "flow" 
associated with the absence of detailed balance becomes 
zero at r — > oo. The effect of the induced sources at the 
special lines and planes can be taken into account per- 
turbatively. The first step is to only consider the effect 
of the origin, as we have done above. The second step 
would be to compute the strength of the sources at the 
lines and planes on the basis of the first order solution, 
and then to determine their function E(k) and recalcu- 
late the solution Eq.(|48|). The induced sources around 
the origin basically lead to a modification of the strength 
M(0), but not of the asymptotic decay. However, the 
far away points at the lines and planes could modify the 
asymptotic decay. A closer inspection of the field of these 
sources shows that it is of order 1/r 7 , since the differing 
rates lead to a local Laplacian acting on the first order 
field decaying as 1 /r 5 . Hence, every step of this perturba- 
tive calculation yields a leading term 1/r 5 ; the amplitude 
changes in every step and its determination is a difficult 
problem indeed. 



D. Correlations for general z 

From the previous section, it is clear that correla- 
tions occur whenever the detailed balance condition is 
frustrated around the origin. The stationary master 
equation then becomes a discrete Poisson's equation in 
(n — l)(d — 1) dimensions, leading to correlations that 
decay with an integer power of the distance r. Following 
the derivation in appendix |b| it is clear that the asymp- 
totic behavior comes from the lowest non-isotropic term 
in -E(k), since division by 1 — D(k) w k 2 gives a sin- 
gularity. The value of the exponent can be calculated 
as 



moments in the fee packing, we find that order 
correlations vanish as 1/r 4 . 



2 and 



(n — l)(d — 1) + order 



(50) 



where (n — l)(d — 1) is the dimensionality of the corre- 
lation space and order is the order of the lowest non- 
isotropic terms in the expansion of E(k). Although this 
result is remarkably simple, the actual calculation of 
E(k) is not trivial, as it reflects the symmetries of the 
jump directions on the (n — l)(d— 1) dimensional lattice. 
Working out the 4-dimensional lattice of the third order 



VI. CORRELATED q-DISTRIBUTIONS 

So far, we have only discussed (/-distributions of the 
type Eq.(^), for which there are no correlations between 
(/-values at different sites. We have shown that, for these 
g-distributions, there are no asymptotic 2-point force 
correlations. In this section we will demonstrate that 
even the smallest correlation between (/-values at differ- 
ent sites induces 2-point force correlations. We first solve 
the problem for arbitrary correlations in the triangular 
packing. Then, we study the fee packing assuming only 
a nearest-neighbor (/-correlation; this already leads to 
force correlations that decay as 1/r 6 . 



A. Triangular packing with arbitrary q— correlations 

In general, the (second order) transition rates are de- 
fined by Eq.(pT|). For z = 2, the displacement vector 
a can only take two values, for which we conveniently 
choose ±|. This allows us to write the transition rates 
as 



W (k) 
W-t{k) 



= % 




-k,+i - 


%- 


-\ Qi+k- 




= % 




-fc,+i 








= % 


2 v 


— Qi+k, 


-|) 


= 1/2- 


W (k), 


= % 




k)— 2 








= % 


r 2 • 


— 1i+k,- 




= 1/2- 


Wo(fe). 



Asymptotically Wo(k) has to approach the value 1/4, for 
(/-distributions without long-ranged correlations. As the 
second moments diffuse on a line, one can easily construct 
a detailed balance solution: 



[1/2 - W (k - 1)] M(k - 1) = [1/2 • 



W (k)]M(k), 

(52) 



= i/2-wb(o; ) 

V ; l/2-W o (A0 V ' 



(53) 



This is the general form of the 2-point force correlations 
M(k) in the triangular packing, as a function of Wo(k) 
that describes the (/-correlations. One can draw two in- 
teresting conclusions from this result. First of all, there 
can only be an uncorrelated solution if Wo(fc) is constant 
(i.e. 1/4) for each k ^ 0. This means that even the small- 
est (/-correlations lead to force correlations. Secondly, the 
long distance behavior of the 2-point force correlations is 
identical to that of the 2-point (/-correlations, following 
from the simplicity of Eq. (J55 
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B. fee packing with nearest neighbor q— correlations 



VII. DISCUSSION 



Unfortunately, the analysis is much more complicated 
for the foe packing, whose second order moments live 
on the 2-dimensional triangular lattice of Fig. |[ We 
therefore allow only correlations between g-values at 
neighboring sites. Remember that one can easily con- 
struct an uncorrelated solution M(r) for uncorrelated 
^-distributions, Eq.([39"|), since all detailed balance con- 
ditions at the origin are identical by symmetry. This 
still holds when there are nearest neighbor correlations. 
However, the detailed balance condition will now be frus- 
trated on the ring of surrounding sites, as these are con- 
nected in four a priori different directions, see Fig. ^|. In 
analogy to the problem discussed in the previous section, 
the stationary master equation for SM(r) transforms into 

5M (r) -1/6^2 5M{r -~f) =p(r). (54) 

The "charge density" p(v) is only non-zero around the 
frustrated ring, see appendix Again, it is a discrete 
version of Poisson's equation, but now in 2-dimcnsions. 
The solution can therefore be expanded in cylindrical 
harmonics, exp(in0)/r", and the six-fold symmetry of 
the lattice requires n > 6. The problem is again solved 



rigorously by Fourier transformation of Eq.(54). In ap- 
pendix O we show that 



SM(r) oc 



cos(6</>) 



(55) 



which is in agreement with the simple electrostatic pic- 
ture. 

o o o 




o o o 

FIG. 5: fee packing: second order moments diffuse on a trian- 
gular lattice. The ring around the origin has differing rates. 



So, for the fee packing, we find that even a nearest- 
neighbor g-correlation leads to 2-point force correlations 
that decay with a power law. This algebraic decay is 
generic for z > 3 since any ^-correlations lead to a mas- 
ter equation whose detailed balance relations cannot be 
solved around the origin. 



We have studied force correlations in the g-model at 
infinite depth, for general ^-distributions. The calcu- 
lated correlation functions are rather unusual: for q— 
distributions of the type Eq.(||), correlations only show 
up at higher orders, and these correlations decay with a 
power of the distance. The only exceptions are the q— 
distributions given by Eq.(|l7|), which do produce a fac- 
torized force distribution. The results for the triangu- 
lar packing and the fee packing are summarized in Ta- 
ble | As an example, consider two different sites i and 
i + k in a layer of the triangular packing. Since there 
are no correlations in the second and third order force 
moments, we find {fifi+k) = 1 and (fff i+ k) = (/ 2 ), 
independent of the distance k. However, the moments 
(fifi+k) and (fffi + k) are correlated and approach their 
asympotitic value as 1/k 5 . The fact that one has to 
go to higher orders to observe force correlations is the 
reason why numerical simulations only marginally differ 
from the mean field solutions [Q . The (single-site) mean 
field solutions p 1 (/) are correct up to the third order 
moments, for the triangular packing. This implies that 
P mt (f) "wiggles" around the real solution p*(f); the de- 
viation p mi (f) — p*{f) changes its sign 4 times (Fig. ||). 

Packings that have more than three g-values per site 
(z > 3) already have third order correlations. Also this 
time correlations only decay algebraically; for the fee 
packing we find 1/r 4 . This algebraic decay can be un- 
derstood from an analogy with electrostatics. The force 
moments evolve according to a master equation, and the 
corresponding stationary state is described by a discrete 
version of Poisson's equation. The "source" turns out 
to be a multipole around the origin, which is only ac- 
tive whenever the master equation has no simple de- 
tailed balance solution. The moments therefore approach 
their asymptotic (uncorrelated) values algebraically. The 
value of the exponent depends on the dimension of the 
correlation space [n — IV <i — I), and on the symmetry of 
the multipole, see Eq.(J50|). 

Although in general correlations do exist, there is a 
special class of g-distributions, given by Eq.([l7]), for 
which there are no force correlations at all. This has 
been demonstrated by means of condition (|l4|), which 
has a nice physical interpretation. It can be shown that 
the function ip(s) is the Laplace transform of the dis- 
tribution of interparticle forces that live on the bonds 
connecting the particles: Vi yCC = qi.afi- Although the g's 
leaving a site are correlated (they have to add up to 1), 
the corresponding Vi, a can become statistically indepen- 
dent. It is only when this miracle happens that the force 
distribution becomes a product state. Nevertheless, the 
q distributions for which this is the case range from in- 
finitely sharp (r — > do) to almost critical (r — > 0). 

Finally, we found that there will be 2-point force cor- 
relations whenever the g-values of different sites are cor- 
related. Even with only nearest neighbor g-correlations, 
the fee packing has force correlations that vanish as 1 / r 6 . 
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Again, the triangular packing is less sensitive for correla- 
tions; the nature of the force correlations is identical to 
that of the g-correlations, Eq. (|53|) . 

Acknowledgements The authors would like to thank 
Wim van Saarloos, Martin van Hecke and Martin Howard 
for stimulating discussions. 



packing 


n = 2 


n — 3 


n = 4 


n — 2, with q-corr. 


triangular 


line 


triangular 


bec 


line 


(d=2) 


no corr. 


no corr. 


1/r 5 


like q-corr. 


fee 


triangular 


4-dirn. 


6- dim. 


triangular 


(d=3) 


no corr. 


1/r 4 




1/r 6 



TABLE I: Summary of the results for the triangular packing 
(z — 2; d — 2) and the fee packing (z = 3; d = 3). The 
nth order force moments diffuse on a (n — l)(d — 1) dimen- 
sional lattice; the lattice structures are listed in the first row. 
The second row shows the nature of the corresponding force 
correlations in the stationary state. 



APPENDIX A: MOMENTS OF 
^-DISTRIBUTIONS 

This appendix is about the moments of the g-distri- 
butions, defined by 



ni no n z 
?1 f 2 



(Al) 

These different moments are not independent because of 
the ^-constraint. As the distributions are normalised, 
the zeroth order moments are unity; the first order mo- 
ments are, of course, all 1/z. All second order moments, 
for which X^ 77 -* = 2, can be described by only one free 
parameter. Defining r\ n as 



In = J r}{q}5 ^1 - qc^j dq 



q t . 



(A2) 



one finds 



X 91® = V2 + (z- 1)<Z1#2 = 

i—l 

J v(q)S - X 9a l dq qi X qi = l l z 



hence 



qi qz 



(z-1) 



(!/*-%)■ 



(A3) 



(A4) 



From a similar argument, one can derive for the third 
order moments 



or = m 



q\qi 



-im - 



(z-l) 

For z — 2 there is even a relation between 773 and r]2 



(A5) 



1 = X q i q o q k = 2 V3 + Gq1q 2 

ijk 



V3= - 1/4. 
For z — 3, there is an additional third moment, namely 



1 = X &?J0ft = 3 % + 18gi92 + 6919293 

919293 = ^(1-9772 + 6773). (A7) 

The extension to higher orders and higher z is straight- 
forward. 

For the special class of 77(g) defined in Eq.([l7|), one 
can calculate the moments r\ n from a generalization of 
Eq.pl) PI 



Vn = 



r(zr)r(r + n) 
r(r)r(zr + 71) 



(A8) 



In order to show that Eq.(|24|) is indeed obeye d by the 
special class (with z = 2), we first invert Eq.(A8) for 
?i = 2 



1 - 2772 

4772 - 1 ' 



(A9) 



From this one can calculate 774 as a function of 772, which 
precisely results in Eq.(p^). A similar inversion for z = 3 
leads to 



1 - 3 ??2 
r = . 

9m - 1 

from which one derives Eq.(|2£ 



(A10) 



APPENDIX B: THE BCC LATTICE 

In the triangular packing, the fourth order force mo- 
ments (fifi+kfi+lfi+m) diffuse on the bec lattice of Fig. 
^, with differing jump rates on special lines and planes. 
In this appendix, we list these rates explicitly and we 
solve the corresponding stationary master equation. 

The jump rates can be calculated from 

W 7 (M,m) = qi,a qt+k,a> qt+l.a" qi+k,a>", (Bl) 

with the z 4 = 16 jump directions 

7= (a' -a, a" -a, a'" -a). (B2) 



As a can take the values ±i, there are two self rates for 
which all cv's are the same. As a consequence, there are 
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14 outgoing directions, namely ±(1,0,0), ±(1,1,1), and 
±(1, 1,0) plus their permutations. The first two are di- 
rections for which three of the four a's are equal, and they 
correspond to the corners of Fig. ||; the third represents 
the ju mps towards the body centers. If all position indices 
in Eq.(Bl) are different, the transition rates are simply 
1/z 4 = 1/16. On the special lines and planes where one 
or more position indices coincide, we encounter differing 
rates. The geometry of the problem is summerized in 
Table || 



from \ to 


(0,0,0) (ft, 0,0) (k,k,0) (k,l,0) (k,l,m) 


origin (0, 0, 0) 

line(c) ( fc <°<°) 
(k,k,k) 

line (b) (k,k,0) 

plane ^ 
(k,k,l) 

bulk (k, I, m) 


A ^ o 2 
1 Q 1 7 I O 1 O 

kit lit iiiiz iiii 2 

(q 2 ) 2 qlWq2 (g 2 ) 2 qlWcn (9T92) 2 

1~ 9 l~9 1~9 1 

4<?I III 111 4<?l92 

111 
16 16 16 



TABLE II: The transition rates W 1 (v) for the fourth order 
master equation. 



with e = 8e c = — 6e&. 

For the large r behavior we need the expansions for 
small k. One finds 

4(k) - 1 - \k 2 + ^ [ fc 4 + 4(k 2 x k 2 y + k 2 y k 2 z + klk 2 x ) } + ■■■ 

(B7) 

and 

E b (k) = 1 - + J_[fc 4 _ (fc^ + fc2 fc 2 + k 2 k 2 )] + . . . 

From these expressions one derives the expansion 
1 - D(k) 24 V 32 



(B8) 



7 pp 1 1.2 JUS 1 U2 Ul 

I K X Ky -T ^y^ Z T S/j; 

8 P 



(B9) 



The first two terms in the expansion are regular and thus 
give rise to short range contributions. The last term 
leads to the asymptotic behavior, by means of the in- 
verse Fourier transform 



dk (klk 2 . + k 2 y k 2 z + k 2 z kl) exp(-ik ■ r) 
Vbz ~k 2 ' 



(B10) 



From this table we deduce the rates e c to the corners 



and £b to the body centers, which occur in relation (46) 
We find 



1 

14' 



£6 



2 2 



1 

14' 



(B3) 



and one easily verifies from the property (71+92 = 1 
that the relation 8e c + 6e& = holds. In general, the 
rates do not obey the detailed balance condition Eq.(^) 
in the elementary loop origin-corner-body center-origin. 
Keeping only the rates in this loop as deviations from 
the bulk leads to equation (|47]). For the definition of the 
two functions E(k) and Dk) we introduce two auxiliary 
functions: one for the contribution of the corners 



E c (k) 



COS -7 h COS -f 



k x -\- ky k z 



k x ky k z 



(B4) 



and one related to the body centers 



1 



E^ik) = — (cos k x + cos k y + cos k z ). (B5) 



This integral can be evaluated by differentiation of the 
well-known 



dk cxpHk ■ r) _ 1 
V B z k 2 ~ Anr 7 1 ' 



where a factor k x in Eq.(BK) corresponds to applying 
d/dx. This leads to expression (BT 



APPENDIX C: g— CORRELATIONS IN THE FCC 
PACKING 



In Eq. (|54|) we formulated the problem for the sec- 
ond moments in the fee packing with nearest-neighbor 
g-correlations. The "charge density" p(r) on the right 
hand side of the equation is the product of the moment 
M(-f), referring to the neighbors of the origin (all are 
the same by symmetry), with a function whose Fourier 
transform is given by 



E(k) 



7'i7 



7-7' 



exp[ik • (7 + 7') 



(CI) 



The ID™ — -y' sire the deviations from the bulk transition 
rates 1 /6. These are only non-zero for the ring of nearest 
neighbors around the origin shown in Fig. |5| 



The two functions D(k) and E(k) are then given as 
D(k) 



E(k) = e[E c (k) - E b (k)}, 



w 

(('2 



W4 = 



-£2, 



Wl 

w 3 



w 5 

£0 - 



2e x + 2e 2 . (C2) 



The equalities reflect the symmetry of the triangular lat- 
(B6) tice. Inserting Eq.(|cT|) into the Fourier transform of 
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Eq.(||) leads to 



u - D(k) 



M(r) =2/3 + M(j)J2~ '-— exp(- ik-v). (OS) 



The consistency equation for M('y) follows by taking r as 
one of the nearest neighbors of the origin. The function 
D(k) is given by 



~ if k x + V3ky k x - y/3k y 
D(k) = - I cos k x + cos + cos - 



(C4) 



and E(k) can be expressed as 



E(k)/6 = e [l - £>(2k)] + 2 £l [l - D'(k)} + 2e 2 [l - D(k)}, 

(C5) 

with the new function 



(C6) 



For the asymptotic behavior of M(r) we must make an 
expansion of i?(k)/(l — D(k)). For the first two terms 
we find 

l-£>(k) V 16 256 



1 k\ Gk^ky + 9k^,ky 



192 



A- 2 



,(C7) 



1 8 128 



1 - D(k) 



1 fcg ~ 6fc^fcg + 9fcgfc4 
288 fc 2 



(C8) 

and the third term is simply a constant. The asymptotic 
behavior is given by Fourier inversion of the first singular 
term in k, i.e. 



dk (k x - 6fc*fc 2 + 9/c 2 fc4) exp(-ik • r) 



V BZ 



k 2 



960cos^) (C9) 



7r r 

This integral can be obtained by differentiation of 
dk exp(— ik ■ r) log(L/r) 



V BZ fc 2 2vr 
where L is the size of the system. 
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